% radon_example
% This m-file was based on one written by Gianni Schena at Univ. Trieste.
% It is a natural place to start because it directly compares matlab's
% radon and iradon routines with an algorithm provided in this suite.

if ~isvar('xmat')
	ig = image_geom('nx', 64, 'ny', 64, 'dx', 1);
	xtrue = phantom(ig.nx)'; % make phantom
	xtrue = max(xtrue,0);
	im clf, im pl 2 4
	im(1, xtrue, 'phantom'), cbar

	% project phantom using matlab way
	na = 180/2;
	angle = [0:(na-1)]/na*180;
	[yi_mat, rad] = radon(xtrue', angle);	% trick: transpose!
	printm('ray spacing = %g,%g', min(diff(rad)), max(diff(rad)))

	sg = sino_geom('par', 'nb', size(yi_mat,1), 'na', size(yi_mat,2), ...
		'dr', rad(2)-rad(1));
	im(2, yi_mat, '"radon" sino'), cbar

	% matlab reconstruction from sinogram
	xmat = iradon(yi_mat, sg.ad, 'linear', 'Ram-Lak', 1, ig.nx)'; % transpose!
	im(3, xmat, '"iradon" recon'), cbar
	im(4, abs(xmat-xtrue), '|error|'), cbar
prompt
end

%
% make system model: use Gtomo2_wtmex if possible,
% otherwise fall back to Gtomo2_strip
%
if ~isvar('G'),	printm 'G'

	if has_mex_jf
		G = Gtomo2_wtmex(sg, ig, 'pairs', {'strip_width', sg.dr'});
	else
		printm 'no mex, so reverting to slow Gtomo2_strip: be patient!'
		G = Gtomo2_strip(sg, ig, 'strip_width', sg.dr);
	end
prompt
end

if ~isvar('yi'), printm 'yi' % make sinogram
	yi = G * xtrue; % forward projector, aka DRR
	im(7, yi, 'sinogram'), cbar
	im(8, abs(yi-yi_mat), 'sino |diff|'), cbar
prompt
end

% reconstruct using em_fbp FBP algorthm
if ~isvar('xfbp1'), printm 'xfbp1'
	fbp_kernel = [1];	% noiseless data, so no filtering
	xfbp1 = em_fbp(sg, ig, max(yi,0), [], [], 'kernel', fbp_kernel);
	im(5, xfbp1, 'em\_fbp'), cbar
	im(6, abs(xfbp1-xtrue), 'error'), cbar
prompt
end

if ~isvar('Rq'), printm 'Rq'
	l2b = 2;	% log_2(beta)
%	C = C2sparse('leak', ig.mask, 8, 0);
%	Rq = Robject(ig.mask, 'beta', 2^l2b);
	Rq = Reg1(ig.mask, 'beta', 2^l2b);

	if 1
		psf = qpwls_psf(G, Rq, 1, ig.mask);
		printm('expected FWHM = %g', fwhm2(psf))
	end
prompt
end

% fix: Caller must do G = G(:,mask(:)) for masked reconstructions.

% CG iterative reconstruction
if ~isvar('xcg'), printm 'xcg'
	wi = ones(size(yi));
	W = diag_sp(wi(:));
	niter = 10;
	xcg = qpwls_pcg(xfbp1(ig.mask), G, W, yi(:), 0, ...
		Rq.C, 1, niter, ig.mask, 1);
	xcg = ig.embed(xcg(:,end));
	im(7, xcg, 'CG'), cbar
	im(8, abs(xcg-xtrue), 'error'), cbar
prompt
end

%
% run the block-iterative SPS-OS algorithm
%
if ~isvar('xsps'), printm 'xsps'
	nblock = 10;
	Gb = Gblock(G, nblock, 0);

	if ~isvar('R'),	printm 'R'
%		R = Robject(ig.mask, 'type_denom', 'matlab', ...
%			'potential', 'huber', 'beta', 2^l2b, 'delta', 0.1);
		R = Reg1(ig.mask, 'type_denom', 'matlab', ...
			'pot_arg', {'hyper3', 0.1}, 'beta', 2^l2b);
	end

	xsps = pwls_sps_os(xfbp1(ig.mask), yi, [], Gb, R, niter, 2, [], [], 1, 1);
	xsps = ig.embed(xsps(:,end));
	im(7, xsps, 'SPS-OS'), cbar
	im(8, abs(xsps-xtrue), 'error'), cbar
end

% compare NRMSE
pr nrms(xfbp1(:), xtrue(:))
pr nrms(xmat(:), xtrue(:))
pr nrms(xsps(:), xtrue(:))
